/*=================================================================================
 * The functions of Euler
 *=================================================================================*/

/*---------------------------------------------------------------------------------
 * Standard C++ library headers
 *---------------------------------------------------------------------------------*/

#include <cstring>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>

/*---------------------------------------------------------------------------------
 * Local headers
 *---------------------------------------------------------------------------------*/

#include "Euler.h"
#include "Global.h"
#include "_string_switch.h"
#include "FileReader.h"
#include "mpi.h"
#include "omp.h"

using namespace std;

/*---------------------------------------------------------------------------------
	Name:		Euler::Euler

	Input(0):	None

	Function:	Read files which identify the physical problem;
				Initialization;
				Compute the flow field and output in real time;
				Output the final result

	Return:		None
 *---------------------------------------------------------------------------------*/

Euler::Euler()
{
	/* Read the ctrl file */
	FileRead();
	
	/*MPI information*/
	MPI_Comm_rank(MPI_COMM_WORLD, &myid);
	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
	num_thread = 4 ;
	if (myid == 0) {
		cout<<"------Beginning of the main part------"<<endl;
		cout << "Input ctrl file is:" << this->ctrlfile << endl;
		cout << "Total process: " << numprocs << endl;
		cout<<"------End of MPI information------"<<endl;
	}
	
	/* Read the mesh file and process the mesh */
	Mymesh.MeshProcess_MPI(myid,numprocs,meshFile);

	/* Read the input file */
	TwoDim.InputRead(inputFile, Mymesh.xnode, Mymesh.ynode, Mymesh.bc);

	/* Read the reaction model and thermo file */
	React.ReactionRead(Reaction_Model, thermoFile, Mymesh.xnode, Mymesh.ynode, Mymesh.bc);

	/* Intitalize the flow field */
	TwoDim.FieldInitial(React.Ri, React.Mw, React.Mi, React.Yi, React.Coeff0, React.Coeff1);

	/*Initialize MPI data */
	TwoDim.MPIdataInit(numprocs);
	double begin_time = MPI_Wtime();

	/* Compute the flow field and output in real time */
	Computing(count);
	
	double end_time = MPI_Wtime ();
	/* Output the final result in a file */
	if (myid == 0) printf("Here is end\n");
	MPI_Barrier(MPI_COMM_WORLD);
	Output_MPI(count,myid);
	if (myid == 0) printf("Total time = %lf", end_time-begin_time);
	return;
}

/*---------------------------------------------------------------------------------
	Name:		Euler::FileRead

	Input(0):	None

	Function:	Read the ctrl file to get the global informartion

	Return:		None
 *---------------------------------------------------------------------------------*/

void Euler::FileRead()
{
	char react[50] = {0}, diff[50] = {0}, timeadv[50] = {0};

	FileReader Ctrl;
	Ctrl.readFile(ctrlfile);

	/* Read number */
	Final_Time = Ctrl.getdoubleParameter("Final-Time");
	cfl = Ctrl.getdoubleParameter("cfl");
	if (cfl == 0)
		dt = Ctrl.getdoubleParameter("dt");
	count = Ctrl.getIntParameter("count");

	/* Read string */
	strcpy(meshFile, Ctrl.getStringParameter("gridname").c_str());
	strcpy(inputFile, Ctrl.getStringParameter("initialization").c_str());
	strcpy(Reaction_Model, Ctrl.getStringParameter("reaction-model").c_str());
	strcpy(thermoFile, Ctrl.getStringParameter("Thermofile").c_str());

	/* Judge the reaction scheme */
	strcpy(react, Ctrl.getStringParameter("reaction").c_str());
	switch (hash_(react))
	{
	case "Trapezoid"_hash:
		Reaction_Sch = Global::Trapezoid;
		type = 0;
		break;
	case "Trapezoid-Adaptive"_hash:
		Reaction_Sch = Global::Trapezoid;
		type = 1;
		break;
	case "IMEX"_hash:
		Reaction_Sch = Global::IMEX;
		break;
	case "Cantera"_hash:
		Reaction_Sch = Global::DNN;
		type = 0;
		break;
	case "DNN"_hash:
		Reaction_Sch = Global::DNN;
		type = 1;
		break;
	default:
		cout << "Invalid reaction scheme" << endl;
	}

	/* Judge the difference scheme */
	strcpy(diff, Ctrl.getStringParameter("difference").c_str());
	switch (hash_(diff))
	{
	case "MUSCL_1"_hash:
		Diff_Sch = Global::MUSCL_1;
		break;
	case "MUSCL_2"_hash:
		Diff_Sch = Global::MUSCL_2;
		break;
	case "WENO"_hash:
		Diff_Sch = Global::WENO;
		break;
	default:
		cout << "Invalid difference scheme" << endl;
	}

	/* Judge the time advance scheme */
	strcpy(timeadv, Ctrl.getStringParameter("time-adv").c_str());
	switch (hash_(timeadv))
	{
	case "EE"_hash:
		TimeAdv_Sch = Global::EE;
		break;
	case "TVD_RK3"_hash:
		TimeAdv_Sch = Global::TVD_RK3;
		break;
	default:
		cout << "Invalid time-advance scheme" << endl;
	}

	return;
}

/*---------------------------------------------------------------------------------
	Name:		Euler::Computing

	Input(1):	The number of output files

	Function:	Solve and update the flow field;
				Output the result in real time

	Return:		None
 *---------------------------------------------------------------------------------*/

void Euler::Computing(int count)
{
	double time = 0.0, TimeX[5]={0.0};
	double WTime_cal[2]={0.0};
	int num = 0;
	iteration = 0;

	while (time < Final_Time)
	{
		WTime_cal[0]=MPI_Wtime();
		/* Calculate the time step */
		if (cfl == 0)
			TwoDim.dt = dt;
		else if (Reaction_Sch == Global::DNN)
			TwoDim.dt = dt / 2;
		else
			TwoDim.CFLcondition_MPI(myid,cfl, Final_Time);
		WTime_cal[1]=MPI_Wtime();
		TimeX[0] += (WTime_cal[1]-WTime_cal[0]);

		/* Ouput the real-time result unifomrly */
		if (time <= num * Final_Time / count && time + TwoDim.dt > num * Final_Time / count)
		{
			Output_MPI(num,myid);
			num = num + 1;
		}

		//TwoDim.SweepX();
		/* Update the boundary */
		WTime_cal[0]=MPI_Wtime();
		TwoDim.FieldBoundary_MPI(myid, numprocs);
		TwoDim.MPIBoundary(myid,numprocs);
		WTime_cal[1]=MPI_Wtime();
		TimeX[1] += (WTime_cal[1]-WTime_cal[0]);

		/* Solve the advection term */
		WTime_cal[0]=MPI_Wtime();
		TwoDim.Advection(TimeAdv_Sch, Diff_Sch);
		WTime_cal[1]=MPI_Wtime();
		TimeX[2] += (WTime_cal[1]-WTime_cal[0]);
		WTime_cal[0]=MPI_Wtime();
		/* Identify the reaction scheme */
		switch (Reaction_Sch)
		{
		case 0:
			Trapezoid(type);
			break;
		case 1:
			IMEX();
			break;
		case 2:
			DNN(type);
			break;
		}
		WTime_cal[1]=MPI_Wtime();
		TimeX[3] += (WTime_cal[1]-WTime_cal[0]);
		TimeX[4] = TimeX[0]+TimeX[1]+TimeX[2]+TimeX[3];
		/* Monitor the iteration */
		time += TwoDim.dt;
		iteration += 1;
		if (myid == 0){
			 cout << time<<"  ";
			 for (int i=0;i<5;i++)
			 cout<<TimeX[i]<<"  ";
			 cout << endl;
		}
		/* Judge if there exists nan value */
		if (TwoDim.D.IsNan())
		{
			cout << "D is nan\n";
			abort();
		}

		if (TwoDim.T.IsNan())
		{
			cout << "T is nan\n";
			abort();
		}
	}
	if (myid == 0) OutputWtime(TimeX);
	return;
}

/*---------------------------------------------------------------------------------
	Name:		Euler::Trapezoid

	Input(1):	The type of the Trapezoid formula

	Function:	Solve the reaction term by trapezoid formula
				(Adaptive (1) or not (0))

	Return:		None
 *---------------------------------------------------------------------------------*/

void Euler::Trapezoid(int type)
{
	/* Update the flow field after solving the advection term */
	TwoDim.Update_after_Adv();

	if (type == 0)
	{
		/* Solve the reaction source term */
		TwoDim.Di = React.Trapezoid(TwoDim.Mc, TwoDim.Di, TwoDim.T, TwoDim.dt);

		/* Update the flow field */
		TwoDim.Explicit();
	}
	else if (type == 1) /* Adaptive */
	{
		for (int i = 0; i < React.NchemNow; i++)
		{
			/* Solve the reaction source term */
			TwoDim.Di = React.Trapezoid(TwoDim.Mc, TwoDim.Di, TwoDim.Yi, TwoDim.T, TwoDim.dt, i);

			/* Update the flow field */
			TwoDim.Explicit(React.Nchem, i);
		}

		/* Monitor the reaction iteration step */
		cout << "Reaction iteration step: " << React.NchemNow << endl;
	}
}

/*---------------------------------------------------------------------------------
	Name:		Euler::IMEX

	Input(0):	None

	Function:	Solve the reaction term by IMEX method

	Return:		None
 *---------------------------------------------------------------------------------*/

void Euler::IMEX()
{
	/* Solve the reaction source term */
	React.Diagonalized(TwoDim.Mc, TwoDim.Di, TwoDim.T, TwoDim.Partial_T);

	/* Update the flow field */
	TwoDim.Update_IMEX(React.CMS, React.MD);
}

/*---------------------------------------------------------------------------------
	Name:		Euler::DNN

	Input(1):	The type of the reaction term solver

	Function:	Solve the reaction term by Cantera or DNN method
				(DNN (1) or Cantera (0))

	Return:		None
 *---------------------------------------------------------------------------------*/

void Euler::DNN(int type)
{
	TwoDim.Update_after_Adv();

	TwoDim.dt = dt;
	TwoDim.Di = React.ReactionS(type, TwoDim.Di, TwoDim.T, TwoDim.P, TwoDim.dt);
	TwoDim.Explicit();
	
	TwoDim.dt = dt/2;
	TwoDim.Advection(TimeAdv_Sch, Diff_Sch);
	TwoDim.Update_after_Adv();

	TwoDim.dt = dt;
}

/*---------------------------------------------------------------------------------
	Name:		NS::Output

	Input(1):	the number of the output file

	Function:	Output the result in a file

	Return:		None
 *---------------------------------------------------------------------------------*/

void Euler::Output(int num)
{
	int ni = Mymesh.ni;
	int nj = Mymesh.nj;
	int bc = Mymesh.bc;

	/* Set the filename */
	string filename = "./output/File00" + to_string(num) + ".dat";

	/* Output the result */
	ofstream outfile;
	outfile.open(filename, ios::out);
	outfile << "Variables = X,Y,T,D,P,U,V,Ma,Y1,Y2,Y3,Y4,Y5,Y6,Y7\n";
	outfile << "ZONE I=" << ni << '\t' << "J=" << nj << '\n';
	outfile << "datapacking=block\n";

	/* The result without boundary */
	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << Mymesh.xnode(i) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << Mymesh.ynode(j) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.T(i, j) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.D(i, j) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.P(i, j) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.U(i, j) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.V(i, j) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Ma(i, j) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 0) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 1) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, bc) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 3) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 4) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 5) << '\n';

	for (int j = bc; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 6) << '\n';

	/* The result with boundary */
	/*outfile << "Variables = X,Y,T,D,P,U,V,Ma,Y1,Y2,Y3,Y4,Y5,Y6,Y7" << '\n';
	outfile << "ZONE I=" << ni+2*bc << '\t' << "J=" << nj+2*bc << '\n';
	outfile << "datapacking=block" << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << Mymesh.xnode(i) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << Mymesh.ynode(j) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.T(i, j) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.D(i, j) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.P(i, j) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.U(i, j) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.V(i, j) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.Ma(i, j) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.Yi(i, j, 0) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.Yi(i, j, 1) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.Yi(i, j, 2) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.Yi(i, j, 3) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.Yi(i, j, 4) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.Yi(i, j, 5) << '\n';

	for (int j = 0; j < nj + 2*bc; j++)
		for (int i = 0; i < ni + 2*bc; i++)
			outfile << TwoDim.Yi(i, j, 6) << '\n';*/

	outfile.close();

	return;
}

/*at process 0*/
void Euler::Output_MPI(int num,int id)
{
	int ni = Mymesh.ni;
	int nj = Mymesh.nj;
	int bc = Mymesh.bc;

	/* Set the filename */
	string filename = "./output/File00" + to_string(num) +"_" + to_string(id) +  ".dat";

	/* Output the result */
	ofstream outfile;
	outfile.open(filename, ios::out);
	outfile << "Variables = X,Y,T,D,P,U,V,Ma,Y1,Y2,Y3,Y4,Y5,Y6,Y7\n";
	outfile << "ZONE I=" << ni << '\t' << "J=" << nj+1 << '\n';
	outfile << "datapacking=block\n";

	/* The result without boundary */
	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << Mymesh.xnode(i) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << Mymesh.ynode(j) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.T(i, j) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.D(i, j) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.P(i, j) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.U(i, j) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.V(i, j) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Ma(i, j) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 0) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 1) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, bc) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 3) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 4) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 5) << '\n';

	for (int j = bc-1; j < nj + bc; j++)
		for (int i = bc; i < ni + bc; i++)
			outfile << TwoDim.Yi(i, j, 6) << '\n';

	outfile.close();

	return;
}

/*---------------------------------------------------------------------------------
	Name:		NS::OutputWTime

	Input(1):	the value of Wall Time

	Function:	Output the result in a file

	Return:		None
 *---------------------------------------------------------------------------------*/

void Euler::OutputWtime(double* val)
{

	/* Set the filename */
	string filename = "./output/WallTime.txt";

	/* Output the result */
	ofstream outfile;
	outfile.open(filename, ios::out);
	outfile << "Time step Time = "<<*val<<"s\n";
	val ++ ;
	outfile << "Boundary Time = "<<*val<<"s\n";
	val ++ ;
	outfile << "Advection Time = "<<*val<<"s\n";
	val ++ ;
	outfile << "React-Update Time = "<<*val<<"s\n";
	val ++ ;
	outfile << "Total Wall Time = "<<*val<<"s\n";
	outfile.close();

	return;
}

